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A quantum theory of dispersion for an inhomogeneous solid is obtained, from a starting point 
of multipolar coupled atoms interacting with an electromagnetic field. The dispersion relations 
obtained are equivalent to the standard classical Sellmeir equations obtained from the Drude-Lorentz 
model. In the homogeneous (plane-wave) case, we obtain the detailed quantum mode structure of 
the coupled polariton fields, and show that the mode expansion in all branches of the dispersion 
relation is completely defined by the refractive index and the group-velocity for the polaritons. We 
demonstrate a straightforward procedure for exactly diagonalizing the Hamiltonian in one, two or 
three-dimensional environments, even in the presence of longitudinal phonon-exciton dispersion, and 
. an arbitrary number of resonant transitions with different frequencies. This is essential, since it is 

necessary to include at least one phonon (I.R.) and one exciton (U.V.) mode, in order to accurately 
represent dispersion in transparent solid media. Our method of diagonalization does not require an 
explicit solution of the dispersion relation, but relies instead on the analytic properties of Cauchy 
contour integrals over all possible mode frequencies. When there is longitudinal phonon dispersion, 
the relevant group- velocity term is modified so that it only includes the purely electromagnetic part 
of the group velocity. 
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I. INTRODUCTION 
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^5 ' Recent developments in quantum optics have led to the realization that it is essential to include dispersion in the 
quantum theory of a linear or nonlinear dielectric, as in a fiber waveguide. At the simplest level, it is clear that all 
dielectric solids have dispersion and absorption. There are an increasing number of experiments that allow direct 
Q ' access to the quantum nature of propagating radiation fields, ranging from quantum soliton experiments in optical 
fibers to QND measurements; as well as quantum dynamics experiments directed at reducing photon number noise 
^ ' levels for broad-band communications, ultra-precise measurements and other device applications Significantly, 
^ , these experiments - as well as more recent theoretical proposals - have the character of fundamental tests of the 
^ ' quantum theory of interacting fields 0- and of quantum measurement theory. These are complementary to older 
Q-i' accelerator-based tests, taking place in low-energy regimes where a considerable precision of control is available on 
IL* the dimensionality of the interacting quantum fields. Measurements that can be used range from the usual particle- 
counting measurements, to interference based techniques that allow an unprecedented level of information about 
quantum phases. Current technologies even allow measurements of electromagnetic properties of Bose condensates. 
^ ■ Because the effects of propagation always involve more than one spatial mode, it is essential to include dispersion 
in any physical model of a linear waveguide in one or more spatial dimensions. Treating dispersion can present 
problems, because it arises from the fact that the response of the medium to the field is not instantaneous. The 
polarization at a given time depends not only on the field at that time but on the values of the field at previous times 
as well. This nonlocal behavior makes standard macroscopic quantization, in which the medium is represented by its 
susceptibilities, problematic. Even without explicit dispersion, some early treatments even resulted in Hamiltonians 
that did not generate the Maxwell equations at all. This problem was analyzed by Hillery and Mlodinow who 
showed that the problem was caused by the use of an incorrect canonical description. 

However, a previous macroscopic model of a nonlinear, dispersive medium resulted in a quantum theory that has a 
consistent Lagrangian . This was accomplished by breaking the field up into narrow frequency-band components 
and quantizing these separately. The frequency dependent linear susceptibilityis expanded up to second order in 
each frequency band and this results in a local Lagrangian in each band. The present paper is intended to treat the 
microscopic basis of the linear dispersion more precisely. Our model is the quantized version of the classical Drude- 
Lorentz [|l3| single-electron oscillator model, in a linearized continuum limit. Similar continuum models have been 
treated previously, generally in the minimal coupling gauge. The present approach uses the multipolar gauge, which 
eliminates complications due to Coulomb interactions and contact (A^) terms. While these effects can be included 
in minimal coupling theories, they result in considerable complexity. The present approach includes all these effects 
by virtue of the transformation to a multipolar interaction, in which the coupling is to the polarization rather than 
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to the electronic momentum. This has the advantage that inhomogeneous media or higher dimensions can be treated 
easily. 

Some earlier dispersive models of a similar type are known, starting from the phonon theory of Born and Huang 
10 and the exciton theory of Hopfield which used minimal coupling. Hopficld considered the electromagnetic 
field coupled to a polarization field which has a single resonant frequency He diagonalizcd the resulting Hamiltonian 
and found that the resulting elementary excitations, polaritons, have a dispersion relation consisting of two branches 
separated by a forbidden frequency band. More recent work has revisited the Hopfield model, though not always with 
all the terms included in Hopfield's treatment jl^Jl^. Other treatments have added losses to the model by coupling 
the polarization field to loss reservoirs [|l8|,^ , still with just one resonant transition. This has allowed the examination 
of the Kramers-Kronig relations in a fully quantized model. Finally, nonlinear generalizations of of Hopfield's model 
have also been studied In these the linear oscillators of Hopfield's medium have been replaced by two- level atoms 
and the Holstein-Primakoff representation is used to develop a systematic expansion of the polarization of the medium 
which includes nonlinear effects. 

These and other studies have elucidated the fundamental cause of the two main types of dispersion. The first, 
excitonic type is due to electronic transitions, and is typically resonant in the ultra-violet. These are due to relatively 
tightly-bound electrons that are localized to atomic sites in insulators, and somewhat delocalized in semiconductors. 
The resulting classical dispersion is rather well described by the harmonically-bound Drude-Lorenz model. The next 
type is due to phonon transitions, mostly in the infra-red. These excitations are less strongly localized, and have their 
own dispersion relation in the Born-Huang model. The dispersion in the transparent region between these absorption 
bands is caused by the off-resonant, virtual transitions of these two main forms of excitation. In effect, a photon 
propagates in a dielectric as a dressed particle, due to the virtual transitions - and resulting polarization - that is 
induced in the solid. To treat this completely, it is essential to include multiple exciton and phonon resonances, in a 
full three-dimensional model. 

The present paper examines this problem using an approach in which the coupling to the dielectric is included 

2^. This has the well-known advantages that the use of atomic sum-rules is not 
required, since the atomic transition probability for off- resonant (high-frequency) transitions is suppressed in this 
gauge P4 I . The new feature presented here is the inclusion of any number of atomic resonances - thus allowing us to 
recover the Sellmeir dispersion equation, which is known to provide an excellent quantitative description of real 
dispersive, transparent media. The use of multiple resonances is essential to the correct description of the transparent 
region with low group velocity dispersion, that occurs between the absorption bands in most known cases. In addition, 
there is a very straightforward procedure for diagonalizing the Hamiltonian - even in the presence of an arbitrary 
number of resonant transitions with different resonant frequencies. 

In the homogeneous (plane-wave) case, we find the expansion of the fields in terms of the quantized modes has an 
exceptionally simple form, in which only an electromagnetic group- velocity correction has to be included relative to the 
usual non-dispersive theory. When there is no longitudinal phonon/exciton dispersion, this correction term is identical 
to that found in previous phenomenological ]l|, p3| , p^ and single-resonance microscopic models |l|Jl5|,0,0@|l|] . As 
an example. Blow, et. al. |23|] based their expansion of the fields on previous work (due to Loudon p^ ) in which it was 
argued, from energy transport considerations, that the group velocity should appear in the expansion of the phonon 
field in terms of quantized modes, when phonon dispersion is present. Blow, et. al. took this result for phonons and 
used it for the case of an electromagnetic field in a dispersive medium. 

Our work here provides a justification for this quantization procedure from a more fundamental point of view 
- a multiple-resonance microscopic model. It is remarkable that the simple inclusion of a group velocity factor is 
sufficient to exactly diagonalize this complex interaction Hamiltonian. The mathematical technique required to prove 
the results involves the use of Cauchy's theorem to carry out the required sums over the different branches of the 
dispersion relation. This is a significant improvement over previous approaches, since in general there is no algebraic 
solution - even in principle - for the polynomial expressions whose roots give the dispersion relations. In addition, we 
find that when there is phonon or exciton dispersion |^^, which is a realistic feature of many dispersive media, this 
procedure is modified in a straightforward way. Instead of the total group velocity, only the relevant electromagnetic 
component is included in the expansion coefficients, suggesting that the diagonalization is intimately related to the 
power transport processes in the solid. The results apply to one, two or three-dimensional environments, although 
only isotropic dielectrics are included here, for simplicity in the derivations. 

II. ONE-DIMENSIONAL MODEL 

We start by considering a straight-forward one-dimensional model, without the complicating features of transverse 
mode structures and longitudinal phonon/exciton dispersion. This simple case illustrates the essential analytic features 
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of our model. It will be generalized to more realistic, higher-dimensional dielectric structures, in later sections. 



A. One-dimensional Lagrangian 



We start by considering a set of dipolc-coupled electronic Lagrangians for oscillators with mass m,^, charge q^,, and 
N (renormalized) resonant discrete frequencies ^^i--- ^vn- These transitions may correspond physically to different 
types of atom, to distinct transitions within each atom, or more generally to some many-body resonance. Thus, 
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{rl - nlrl) + —q,r,D(x,) 



(2.1) 



Here is the displacement of a charged particle (typically, an electron) in a multi-polar Drude-Lorentz model, from 
the corresponding center of charge (typically, nuclear) position Xi,. Generally, these are all distinct atomic resonances. 
For simplicity, the self-energy terms proportional to , are incorporated into the definition of the resonant frequency 
fl^. The coupling in multipolar gauge is to the displacement field Dix^) at the central position. All Coulomb terms 
in this gauge are carried by the photon-exchange process ||2l[ . 

Next, in order to develop a simple electromagnetic Lagrangian in the multipolar gauge, we introduce a gauge field 
- the dual potential A, so that D = dA/dx and B — ndA/dt. This simply corresponds to a canonical transformation 
to the multipolar gauge ^0 21| of the more usual minimal coupling theory - and is simplified here by the assumption 
that there are no free charges. The discrete atomic positions are replaced by a continuum polarization field, together 
with the appropriate Lorentz-shifts in the resonant frequencies v, due to local field corrections. This simplification is 
used here because it is not really necessary for us to consider the details of local-field corrections at this stage. We will 
show that this approach is able to generate the identical (low-frequency) Sellmeir dispersion equation that is obtained 
in the usual Drude-Lorentz theory, which does include local-field corrections. Of course, the approximations involved 
are only valid for carrier wavelengths much greater than the inter-atomic spacing. 

The corresponding Lagrangian density that generates the correct electromagnetic energy and Maxwell's equations 
in one dimension is, for a one-dimensional wave-guide with effective cross-section A: 



C/A^^AHx) + - 
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P{x)d^A{x) - i(a,A(x))2 



(2.2) 



Here Pi,{x) is the density of the oscillators with resonant frequency and r^{x) is regarded as a continuum field, 
with polarization P{x) — "^yPvix) — '^^qi>Pi>{x)r^{x), in order to allow us to use a continuum approximation in 
later calculations. 

The coupling between the field and the polarization is entirely included in one term (x) , which effectively combines 
all the relevant information about charge, density, and dipole-moment: 



di^ix) = qlp„{x)/{m^eo) 



(2.3) 



Noting that D{x) — eQE{x) + P{x), the electric field is given by E{x) — [D{x) — P{x)]/eo, and the resulting Maxwell- 
Bloch equations have the expected form of: 

d^A - c^dlA = -c^d.^P{x) 

+ ^Iv. = a.dxA . (2.4) 

In order to quantize the system, we simply introduce commutators for the canonical momenta and position, where: 



H(x) = fiA{x) = B{x) 
Tr^{x) = pi,{x)/{eag^{x)) 



(2.5) 



The quantization rules are then the usual ones, obtained on replacing Poisson brackets with commutators, except 
with the atomic operators treated as a continuum field. All position-type operators of the same type of variable must 
commute amongst themselves, as do all momentum-type operators. We list the commutators involving cross-terms 
between the position and momentum operators, and between the different types of variable, for reference in the next 
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sections. In order that the commutators will have a familiar appearance, they are written in terms of the electric 
displacement and magnetic fields, rather than the canonical variables and their momenta. Since the electric field only 
differs from the displacement field by the polarization terms (which commute with field operators at equal times) there 
is no essential difference between the electromagnetic commutators written with the displacement field or the electric 
field. The fact that the electric displacement is the derivative of a potential is, of course, the origin of the differentiated 
delta function in the commutator between the electric and magnetic fields, which is also found in minimal coupling 
theory: 



(I) 


D{x) , B{x') 


= ihS 


(II) [p,{x),^,,{x') 


— ihS 


(III) 


D{x),%{x') 


= 0, 


(IV) 


B{x),n,{x') 


= 0, 


(V) 


D{x),Tf^{x') 


= 0, 


(VI) 


B{x),%{x') 


= 0, 



(2.6) 

it is important to notice here that the commutators are essentially identical to those for the corresponding free fields 
^ , p6| and oscillators. This is a consequence of the fact that the couplings do not involve time-derivatives. It also 
demonstrates that the present theory is canonically equivalent to other techniques involving minimal coupling pO| , pl| . 
One apparent difference is in the commutator between the displacement field and the momentum jf^, , which replaces a 
commutator between the electric field and a gauge-dependent canonical momentum in minimal-coupling theory. Since 
the electric field and displacement field differ by a term that doesn't commute with the canonical momentum, it might 
seem that this introduces a difference. The explanation is due to the well-known fact that the canonical momentum 
in this case is just the usual mechanical momentum, and therefore differs from the minimal coupling momentum - 
which includes a term proportional to the vector potential A. These two effects cancel each other, so the overall 
commutators are unchanged. 

The final Hamiltonian density, equal to the linear energy density of the coupled system, has a rather straightforward 
expression in which there are no explicit Coulomb interaction terms. This is a typical property of multipolar interaction 
Hamiltonians. Effectively, the Coulomb couplings are included partly in the oscillator self- frequency terms (i.e., in 
Vi), and partly in the photon-exchange dynamics that are implicit in the Hamiltonian. The result is that: 



n/A = ^D\x) + ^B\x) - ^D{x)P{x) 
2£o 2/^ eo 

1 ^ 
2^09 Ax) 



(2.7) 
(2.8) 



All these equations hold for an arbitrary spatial distribution p{x) of the polarizable atoms. If p(x) is a sum of 
delta-function terms, a discrete lattice model is obtained; it is an unusual waveguide in which the atoms are all 
lined up in a one-dimensional row, but not an impossible one, with atom-trapping techniques. For simplicity, we do 
not treat this type of model here. Instead, we will focus on the uniform dielectric limit, in which all the local-field 
corrections are included via the relevant Lorentz shifts of the oscillator frequencies, to give a more tractable case. 



B. Mode Structure 



We now specialize to the case of a continuum model with a uniform distribution, in order to find the dispersion 
relations for plane-waves. We introduce a mode structure by finding the eigenmodes of the equation of motion. 
Suppose the solutions to Maxwell's equations have the form: 



X{t,x) 
p^{t,x) 



A it,x) - 

It follows that these satisfy the resulting equations in the form 



A \ ^ikx—iujt 



(2.9) 
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Solving for p^, we find that: 



(2.10) 



(2.11) 



or, summing over all the oscillators and introducing p — 'YliyPv, we find: 
Eliminating the polarization field p leads to the eigenvalue equation: 



(2.12) 
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(2.13) 



We find a band-structure in which there are typically (A^ + 1) positive roots Lu^{k) with /i = 0, 1, ..N to the solution. 
To look at this differently, we can solve for the wave-number fc, at any given frequency uj: 



k = ± 



1/2 



±k{uj) 



(2.14) 



This solution is unique for any given modal frequency, but has forbidden regions at w ~ 0,,^, where fc^ — > —oo. This 
indicates a resonance, or absorption band. In the transmission bands, there is a unique refractive index n{u), and 
hence a well-defined permittivity, e{uj) = n'^{Lo)eQ. 

It should be realized that the dispersion relation is not completely identical to the usual classical Sellmeir expansion, 
although it is very similar to it. The Sellmeir expansion is: 



n^{uj) = 1 + 



9n 



(2.15) 



This agrees with the functional form of the present result to lowest order in g^^ if we assume that gn — gp,. The 
difference is simply due to different approaches to treating local-field corrections in the continuum approximation. If 
a precise local-field theory is required, then the assumption of a homogeneous polarization field must be replaced by a 
lattice model. In fact, the usual Drude-Lorentz model is not typically obtained from a coupled Lagrangian, so it cannot 
be readily quantized directly. Rather, it is obtained from an approximate theory in which the local (microscopic) 
E field plays the role that the D field does in the present theory ||l^. From microscopic considerations, both forms 
needs to have local-field corrections included near an absorption band, in order to give an accurate comparison with 
a three-dimensional crystal lattice, from first principles. 

When this is done in the Drude-Lorentz model, all the low-frequency resonances are shifted by amounts known as 
the Lorentz shift. With these shifts included, the Sellmeir expansion is obtained, with local-field corrections included. 
However, the number of poles in the rational function representation derived here is finite, just as in the classical 
Drude-Lorentz theory. Accordingly, it is always possible to re-express our dispersion relation exactly in the Sellmeir 
rational-function form, using partial fraction expansions, just as in the Drude-Lorentz derivation of the Sellmeir 
equation. In this case, the frequencies and couplings appearing in the final Sellmeir formula are not identical with 
the original frequencies in the Lagrangian; these shifts of course, can be regarded as evidence of the photonic coupling 
between the oscillators in our model. 

An important, and experimentally well-tested property of the Sellmeir equation is that the refractive index ap- 
proaches unity at high enough frequencies, while at low frequencies it approaches a constant value greater than 



lim n {uj) 



(2.16) 
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Our dispersion relation from the multipolar Hamiltonian has a similar behavior, except that the algebraic form is 
slightly different at low frequencies: 

-1 

(2.17) 

Clearly, one form can always be transformed into another, given the obvious restrictions on having distinct roots, 
with X^i/ fi'/^^ ^ 1- We note that this equivalence does not apply in all other models of the dispersion relations, 
which may have inequivalent analytic properties - leading to a different limiting behavior at high and low frequencies. 
Provided the general analytic properties are equivalent, the partial- fraction procedure to transform one form into the 
other is not required in most cases. We shall demonstrate that only the refractive index and group velocity are needed 
to obtain a complete quantum theoretic description of the modes. This information is readily available from the usual 
Sellmeir parameters that are experimentally known for many transparent materials. A typical dispersion relation for 
a solid with both low and high-frequency resonances, would show multiple transmission and absorption bands - with 
three distinct branches to the dispersion curve - and a region of relatively low group-velocity dispersion between the 
absorption bands. This is the origin of the well-known zero-dispersion point (at around A = 1.5/xm in fused silica), 
which plays a prominent role in optical communications systems. 



hm n (uj) 
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III. ONE-DIMENSIONAL MODE OPERATORS 



Having derived the modal solutions, we now wish to expand the fields of the theory in terms of annihilation and 
creation operators. We develop this expansion in two stages. First, we will consider the conditions on the mode 
expansion which preserve the correct commutation relations for the original canonical fields. Then, we show that a 
mode expansion defined this way does lead to a diagonal Hamiltonian form, when the Hamiltonian is re-expressed in 
terms of the mode operators. That is, our goal is to find operators a^{k) which have the function of diagonalizing the 
Hamiltonian, giving the final structure of: 

N 

H = Y, niVf,{k)al{t, k)a^{t, k) dk . (3.1) 



A. Commutation properties 



Clearly, since the above expansion must lead to the same eigenfrequencies as the original Maxwell equations, we 
should define mode operators relative to each branch of the dispersion relation (with, for the sake of simplicity, 
A'' -|- 1 distinct branches) so that: 



A{t,x) = ^ dk [A^(A;)a^(t, fc)e*'=^ + h.c] 



(3.2) 



Here w^(fc) is the inverse of k{u)), for the /z-th branch. The summation is over the discrete branches in the dispersion 
relation, each of which correspond to a different 'particle' type. The time-dependence of the mode operators in the 
Heisenberg picture - given the desired Hamiltonian structure - must be: 



an{t,k) =aij,{k)e 



(3.3) 



These combined field-polarization excitations can be termed polaritons. We will suppose that the a^lk) are chosen 
so that: 



a^{k),al,{k') =5^^,5{k-k') . 



Similarly, the momentum field can be expanded as: 



\l{t,x) = J2 dk n^(/e)a^(fc)e*'='^-*'^''W* -h/i.c. 



(3.4) 



(3.5) 
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The requirement of commutation relations means that we must have (at equal times): 



A(a;),n(a;') =inS{x-x')/A = ^ dk A,,(fc)n;(fc)e*'=(^-^ ^ 
This implies the relationship that, in order to preserve commutation relations, 

N 



h.c. 



(3.6) 



(3.7) 



The Lagrangian density for this model implies that 11 = fidtA. With the assumed time dependence of the annihi- 
lation operators, we can also write n^(fc) in the form of: 

n^(fc) = ~iujf,{k)fiA,,{k) . (3.8) 

The above result therefore reduces to an equation for the expansion coefficients Ap(fc), in the form of: 

N 



Y^uo,{k)Kl{k)^ 



(3.9) 



Next, we wish to obtain an expression for A^i(fc), which is unknown at this stage. It is no restriction to choose 
Ap(fc) to be real. Therefore, we can always choose an (unknown) function w^(fc) so that, in analogy to the standard 
vacuum expansion, 



^Vf,{k)e^{k) 



inAk 



1/2 



(3.10) 



If Vfj_{k) = c and s^{k) = Eq, this reduces to the accepted vacuum field expansion. More generally, we define 
e^(fc) = fc^/[/ia;^(fc)] as the effective permittivity of the /x-tli branch. We will show later that w^i(fc) must be interpreted 
as the group velocity, just as in an earlier narrow-band analysis of this problem, using effective Lagrangian arguments 

In order to demonstrate this, we first substitute the above expression for A^(fe) into the equation for the consistency 
of the field and mode-operator commutation relations (i.e. for [A, II] and [a, a''']). This leads to the very simple result 
that: 



JV 



(I) 



LO^{k) 



(3.11) 



As explained above, we have to determine a function w^(A;) which satisfies this condition, and we intend to demon- 
strate that the choice of v^{k) as the group- velocity of the relevant polariton branch is sufficient; no other correction 
factors are needed in this simple model. In order to verify this, we can differentiate both sides of Eq. ( 2.13| ) with 
respect to k. This gives a group velocity of: 



duj^jk) _ LUf^{k) 
dk k 



1+E 



(3.12) 



which is the function we propose to use in the mode expansion. 

However, it is clear that the mode function expansion of and tTi, are also needed, for a complete demonstration 
of consistency. Using Maxwell's equations, we note that for a Fourier component of at frequency oj, wave- vector 
fc, we must have: 



ikg^X 



(3.13) 



Thus, if we expand Pi, as: 
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N 

E 

/i=0 



dk 



pl{k)a^{k)e 



ikx—ioj^{k)t 



+ h.c. 



it follows that the expansion coefficient for in the j-th frequency band is: 

^ ifcg.A^(fc) 
P'^^"' (f22-a;2(fc)) • 

Similarly, if the canonical momentum for the atomic polarization field is: 

N ^ 

f<(fc)a^(A;)e*'=^-'"''('=)* + h.c. 



-^1 



dk 



then the corresponding expansion coefficient is: 



fcu;^(fc)A^(fc) 
eoini - ul{k)) 



(3.14) 



(3.15) 



(3.16) 



(3.17) 



For these operators to have the correct equal-time commutators, the different oscillator position operators must 
commute amongst themselves at equal times, as must the different momentum operators 7r,y. This is trivial from the 
form of the mode operator expansion. However, the commutation relations (II) between the position and momentum 
operators are non-trivial. The relevant commutation conditions are therefore: 



\$^{x),T:^>{x')]=iM^^>5{x-x')/A = E / 

n 

This in turn implies that: 

N 



dk 



p-;{k)nf{k)e 



ik{x—x') 



h.c. 



(3.18) 



N 



ih 

AttA 



(3.19) 



Expanding the coefficients gives the new equation: 

^ uJ,{k)k^9.Kl{k) 

earn 

fj.= 



However, since £^(fc) = c fc £o/[w^(fc)], and hence 



this simplifies to the form: 



(11) 



2.,^ _ nc^epky^jk) 
V^"^" AwAuf,{k) ' 

c^k^Vt,{k)gv 

u^^{k){u;l{k) - 02)(a;2(fc) - 02,) 



JV 

E 



(3.20) 



(3.21) 



(3.22) 



Finally, to ensure that there are correct field- atom commutators, we must satisfy the commutators HI- VI. For these 
cross-terms between the oscillators and field variables, we find that conditions (III) and (IV), involving commutators 
between the field and the particle position (or the field momentum and particle momentum) are automatically satisfied. 
This occurs for the same reason that commutators like [A(a;), A(a;')] or [II(a;), n(a;')] must equal zero in our mode 
expansion. In all these cases involving pairs of canonical position-type operators or pairs of momentum-type operators, 
the commutator reduces to an integral over fc, which is an odd function oi k - and hence vanishes when integrated 
over all positive and negative k-valucs. 

This leaves the requirements (V) and (VI) , which are that A and liu vaxisi commiite at equal times , as well as H 
and 6,^. These two requirements both imply the same restriction on the expansion coefficients, and hence on v^{k), 
which is that for all k and v we must have the conditions: 
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Despite the complex nature of each of these conditions - which involve sums over all the roots of the dispersion 
equation, and must be satisfied for all the resonant frequencies v, as well all momenta k - we will show that each of 
these sums can be analytically evaluated without requiring an algebraic solution for the roots. 

B. Conditions on expansion coefficients 

From the previous results, we have shown that the condition on the expansion coeSicient of mode operators is that 
we must find a function v^ik), such that: 

together with an orthogonality condition: 

In addition, to ensure commutation between the particle and electromagnetic fields, we should impose the condition: 



(1. 



__^ul{k){ul{k)-ill) 



Earlier work |l2[ on more phenomenological narrow-band quantum models of dispersion led to the conclusion that, 
for an expansion of modes to be consistent with both Maxwell's equations and the known dispersive energy, it is 
necessary to choose w^(fc) equal to the group velocity. Thus, the use of v^{k) = duj^j,{k)/dk is an obvious choice, but 
it is necessary to demonstrate that this still leads to a complete set of consistent commutation relations. 

C. Single Oscillator Case 

As an example, we consider the single-oscillator case, where the dispersion relation can be treated algebraically as 
the solution of a quadratic equation. In this case the refractive index is given by: 

nM^ = [l-^]-. (3.27) 

In order to show how the Sellmeir form can be regained, we define a new frequency (l^ — il^—g. As long as flf, > g, 
the above equation is equivalent to a Sellmeir type of dispersion relation, which is simply: 

Either equation leads to a quadratic for uj'^, having the form: 

u'^ -Lo^ic^k^ + + nl - g ^0 . (3.29) 
The resonant frequencies at any given wavenumber fc, are then: 



2 

" 2 



We now identify Vf^{k) with the group- velocity on each of the two branches, by taking derivatives with respect to k. 
Thus, assuming g (to have distinct branches): 
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kv±(k) l(c^k^\[ ^P-nl + 2g\ 

±^ ^ ' ^ ll± ^ , (3.31) 



where the quantity A is given by 



A = V(c2fc2 + - 4c2fc2(f}2 _ (3 32) 

Clearly it is necessary to have $72 > in order to have distinct real branches to the dispersion relation, each with 
positive frequency uj. This is precisely the condition under which the Sellmeir expansion is valid, as one might have 
expected. 

Summing over the two branches, we note that (defining K = c^k'^ + ill): 

^ kv±{k) K + 2ig-nl) \ KtA 

On taking the sum, this reduces to the required result of: 

^<^^-E^ = l- (3-34) 

This indicates that the use of group-velocity expansion coefficients appears correct in this case, although we have 
only calculated the simplest of the commutators. However, this algebraic technique is rather clumsy to use for the 
other identities. Even worse, it is not able to be used at all for an arbitrary number of branches; in these more general 
cases there is no closed form expression for the solution to the dispersion equation. 

IV. ANALYTIC PROPERTIES OF DISPERSION RELATIONS 

For the other, more complex, commutation relation identities - or more oscillators - it is preferable to use techniques 
from complex function theory, which transform the sums over roots of the dispersion relation to complex contour 
integrals of related meromorphic functions. The dispersion relations considered here have the general structure of: 



u;2 



This can be written in the form of roots of a polynomial in z — oj^, so that: 



l-Eo/^- (4.1) 



Ka(z^) = h{z^,) , (4.2) 



where = cj2^ k — c^k'^ and: 



aiz) = n(^' - ^) - E ff-^' n - ^) ' (4-3) 
h(,z) = zYi[^l-z) . (4.4) 

V 

Next, in order to test the assumption that the expansion coefficients involve group velocities, we must consider the 
slope of the dispersion relations: 



Hence, 



Ka'(z) + ^a(z) = 6'(z) . (4.5) 
oz 



_ du;^{k) c^kajz) 

"^^^ dk uj^{k){h'{z)~Ka') ■ ^ ' 



A simpler way to write this expression - entirely equivalent to the above definition - is Eq. (3.12). While this form is 



more transparent, the above expression is a rational function, which is amenable to analysis using Cauchy's theorem. 
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A. Condition I 



The sum S^^^ has the form: 



c(I) ^ V = V 



Ka{z) 



z{b'{z)-Ka'{z)) 



Next, consider the complex function, defined for the analytic continuation of z into complex values: 

Ka{z) 



z [6(2;) — Ka{z)] 



(4.7) 



(4.8) 



This generally has {N + 2) poles; and has the property that lim f'^^''{z) ~ l/z"^. For example, in the trivial case of 

|z|^cx> 

no oscillators (A^ = 0), we find that: 



z{z — k) 



(4.9) 



In this case, the identity (I) is satisfied trivially, since it reduces to Res[/^^)(z = c^/c^)] = 1. The sum of residues of 
f^^^z) is zero in this case, which must be true in general for a meromorphic function behaving as f^^\z) ~ 1/z'^ as 
\z\ 00. As usual in complex function theory of the inverse variable (u = l/z), a contour integral around all the 
finite poles turns into a contour integral around zero poles at infinity, and hence must equal zero. Thus, we have the 
general result that: 



= ^Res /W(^) =-l+E 



N+l 



Ka{zij) 



Zfi ib'{z^) - Ka'{z^)) 



(4.10) 



However, this is precisely condition (I), for the A^-oscillator case, since 

N+l 



N+l r 

ft 



Zn [b'izij,) - Ka'{zij,)) 



E 



oj^t{k) 



1 . 



(4.11) 



Thus, the use of complex function theory shows that (I) is always satisfied, provided there are {N + 1) distinct roots. 



B. Condition II 

Similarly, we can prove the other relations. For example, to prove relation (II) we consider: 



f^iz) 



iz-ni)iz-nl,) ■ 

Summing the residues of this function, and noting that lim f^^\z) = —l/z, we find (for v ^ v'): 



(4.12) 



Examining the RHS of the required identity, we must obtain the value of f^^\fl^), evaluated at each resonance: 



(4.13) 



(4.14) 



However, b{fll) = at each resonance, so that f^^\fl1) = — l/fi^. Hence, the RHS of the required identity becomes: 



+ 



+ 



. 



(4.15) 
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In the case that = i^', a double pole is found, so the residue is obtained on differentiating /"^'^(z). We can perform 
this operation most simply in the neighborhood of the double root at 2 = Q^, by using the definition of f''^\z) to 
show that: 



1 



Thus, as ^; — > O^, we find the second term is dominated by the pole in the denominator: 

lim /(i)(z) = -l + ^t£ . 
z^ni z Kg„ 



Hence, on differentiating to obtain the residue, 



1 1 



This is sufficient to complete the proof of the second relation, which is 



(4.16) 



(4.17) 



(4.18) 



(4.19) 



C. Condition III- VI 



As shown previously, the conditions (III) -(IV) are straightforward consequences of the general type of expansion 
chosen here, while conditions (V) and (VI) reduce to an identical summation identity. To obtain this last identity, we 
can now introduce a third analytic function. 



/(I) 



(4.20) 



As well as the poles at 2; = and the (A^ + 1) roots of the dispersion relation, this has a pole a.t z — fi^. However, 
the residues at z = and z = cancel each other, so the sum over the remaining zeros must give zero, as required. 
In summary, we find that summing over the residues gives: 



S^y=J2Res[fi'''\z)]=0 
This proves the last sum-rule requirement on the commutators. 



(4.21) 



V. HAMILTONIAN 

We now wish to show that when the Hamiltonian is expressed in terms of the operators a^(fc) and aj^{k), jjL = 
0, . . . , A/', it is of diagonal form. Our first step is to derive an orthogonality relation for the classical modes. This will 

allow us to show that the Hamiltonian contains only terms of the form a|^(fc)ap(fc). The next step is to examine the 
normalization of the modes. Once this has been done we find that the Hamiltonian is given by 

N 

H = J2 dknuj^{k)dl{k)a^{k). (5.1) 
In order to find the proper orthogonality relation for the modes, let us first define the A' + 1 component vector: 

or Ao = A and Ag = p^^ for s > 1. The eigenvalue equations can be expressed in the form (for each value of k) 
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MX ^uj^X, 



(5.3) 



where the {N + 1) x {N + 1) matrix M is given by 



/ fc^c^ ikc^ ikc^ • ■ ■ \ 



M 



-ikg^^ nl 



V 



(5.4) 



The matrix M is not hermitian, but if it is multipUed by the positive, diagonal, {N + 1) x [N + 1) matrix G, 

Grs = SrsGss where Gss = S ^2^^^ ^ s > 1 ' ^^'^"^ 



then the combination GM is hermitian. Therefore, if 



MX = uJiX 



MA ~ UJ2X , 



then 



This imphes that if uif 7^ , then 



(gma^^^|a<^^)^..,^(a(^Va*'V 



{x''\gx^'^) = o, 



(5.6) 

(5.7) 
(5.8) 

(5.9) 



and we have the desired orthogonahty relation. Expressing this in slightly more generality, we note that for each value 

of k there are + 1 eigenvectors A^^'', /i = 0, . . . , iV, each corresponding to a different eigenvalue Wp(fc). As a result 
we have that 



(5.10) 



We now express the fields in terms of the eigenvectors and substitute them into the Hamiltonian, which is given by 
integrating the Hamiltonian density in Eq. (|2.8D over x. In particular, we have that 



Am = Ao 



n, 



-i[iuj^X\ 



This allows us to use Eqs. (^|^) and ( |3.10 ) when calculating the Hamiltonian, and we find that 

N „ 

H^A^IiAY^ / dk(X^'"\k)\GX^^\k))Lol{k)al{k)a^{k). 



(5.11) 
(5.12) 

(5.13) 



In order to show that the Hamiltonian assumes the form given in Eq. (5.1), and to justify our assumption that 
ia^ — uj^a^, we need to prove that 



4TrfiA{x''^\k)\Gx''^\k))ujf,{k) = h. 



Noting that 



X^^' - A 



/ f^v^lk 



(5.14) 



(5.15) 



we see that Eq.(5.14) will be true if 
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E 



kv,, 



,2\2 



= 1. 



(5.16) 



This implies t hat that v^^ = duf^/dk, from Eq. ( |3.12| ). We can then conclude that the expression for the Hamiltonian 

given in Eq. (5.1) is correct. 

Summarizing, our theory of a linear medium with N resonances is described by the Hamiltonian in Eq. (5.1) and 
the corresponding field operators have the expansion: 



D{t,x) 



E{t,x) 



dk k 



dk 



hs{uj^{k))v^{k) 



AnkA 
hkv^{k) 



1/2 



47re(w^(fc))A 



1/2 



(5.17) 



(5.18) 



B{t,x) 



dk 



47rA 



1/2 



(5.19) 



Here we have also included the electric field expansion ( which is obtained by including the polarization term), 
for comparison with more familiar results. As one might expect, the only difference between the electric field and 
displacement field expansions, is a factor proportional to the dielectric permittivity e{Lu^(k) in each branch of the 
dispersion relation. It is important to notice here that the two fields cannot be related by one, frequency-independent 
coefficient. This is a natural consequence of dispersion, and also occurs in the corresponding classical theory. 

In order to provide a more physical understanding of this result, the summation over the branches in one dimension 
can be replaced by an integral over the mode frequency, in each propagation direction, i.e., define: 



D{t,x) = &+'>it,x)+&-\t,x). 



(5.20) 



Now, since the mode frequency has a bounded range for each root, we can define a frequency dependent mode operator 
as: 



l{uj) = af,{k)/^\Vf,{uj)\ 



(5.21) 



where the appropriate root fi is chosen in each case to correspond to the mode frequency argument - except, of course, 
in the forbidden bands. The commutators of the new mode operators are: 

[a{uj),a^uj')] = 6iuj - oj'), (5.22) 

and the mode expansion is now the same as it would be in a non-dispersive case - except that no modes exist in the 
forbidden bands: 



poof 

5(±)(i,a;) = ±j / duj 
Jo 



fik{uj)e{uj) 



AnA 



1/2 



5(w)e±^'=('^)^-*"* + h.c 



Similar equations hold for the other fields; for example, the electric and magnetic field expansions are just: 



pQCf 

&^\t,x) dio 
Jo 



fik{ijj) 



-,1/2 



and: 



B^^\t,x) 



dio 



4:ne{uj)A 



a{uj)e 



a{uj)e 



±ik{iLj)x — iujt 



h.c. 



(5.23) 



(5.24) 



(5.25) 



A'Ke{uj)A 

. The important point is that wc can exactly absorb the group velocity factor into the frequency integral - which, 
however, is only defined in the range of allowed mode frequencies. This result (also obtained in earlier narrow-band 
Lagrangian approach [l|Jl^, and in a single-resonance model [l^ ), was most clearly emphasized in the single- resonance 
theory of Huttner and Barnett ||l^]. It implies that the two-time correlation function for narrow-band fields in the 
transmission band, is essentially identical to those of the corresponding vacuum fields, apart from the usual reflectivity 
factors. This is a necessary ingredient of any theory of the interface properties of the fields, and will be explored 
in more detail in a subsequent paper. The above mode expansion has been widely used in quantum optics, and the 
present result shows that it is exact for a realistic, multiple-resonance model of a dispersive medium - provided we 
recognize that there are no modes in the forbidden bands. 
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VI. HIGHER-DIMENSIONAL MODELS 



We next consider an n-dimensional Lagrangian for oscillators with mass m^, displacement r^^ effective charge q,j, 
and oscillation frequencies about their center of charge position : 



Z Eq 



(6.1) 



Here [rj^g^] is the dipole moment of a charged particle in a multi-polar Drude-Lorentz model. In this general case, the 
labels i> may correspond either to distinct resonances of one atom or to distinct positions. Each resonance has its own 
corresponding effective charge - and hence dipole moment. Any sum rules are incorporated into the definitions of the 
effective charges involved in a given transition. For simplicity, the self-energy terms proportional to are included 
in the definition of the resonant frequencies, which are defined to diagonalize the individual charge-cell Hamiltonians 
in the multipolar gauge. The coupling in multipolar gauge is to the displacement field D(x^) at the central position 
x,y, used as an origin for defining a local polarization. All inter-atomic Coulomb terms in this gauge, are carried by 
the photon-exchange process, which couples atoms at distinct spatial positions. 

We introduce a vector gauge field - the dual potential A, so that D = V x A and B = fidA/dt. The discrete 
cell positions are now replaced by a continuum polarization field as before, together with the appropriate local field 
corrections. To account for more general dielectric structures that may have local interactions not included in the 
Coulomb corrections, we now include a quadratic dispersion term to describe the residual phonon and exciton 
dispersion [ p5[ , that exists in the absence of long wavelength electromagnetic couplings. 

The simplest Lagrangian density that generates the correct electromagnetic energy and Maxwell's equations for an 
n-dimensional wave-guide with effective cross-section A c?'^~" is: 



2fi eo 



P(x) ■ V X A(x) - i(V X A(x))2 



(6.2) 

[P' W - nlpU^) - a,(x)[V,p,(x)]2] /(2£o5.(x)) . 



N 



Here the polarization density due to all the dipoles is P(x) = Piy(x) — r^qi,pi,{'x.), where Pi/(x) is the number 
density of the oscillators with resonant frequency il^. The partial polarization Pi/(x) is regarded as a continuum field, 
with = 1, ..TV labeling the bare frequency of the elementary phonon and exciton resonances, in the absence of 
coupling to the long wavelength photons. We have assumed that the dispersion of phonon and exciton modes are the 
same for longitudinal and transverse modes. The transverse dispersion will ultimately be modified by the coupling 



between the field and the polarization, which is entirely included in the term g^(x), as defined in Eq. (2.3). 

We impose the usual gauge constraint that V • A = 0, so that the field variable only has transverse degrees of 
freedom; this does not apply to the polarization. In addition, we can impose wave-guiding conditions that A is 
restricted to a one, two or three dimensional manifold. In practise, dispersion occurs in the electromagnetic mode 
functions (which are wavelength dependent) , so that it is necessary to solve for the complete three-dimensional mode 
structure to rigorously treat a fiber waveguide, for example. However, a simple low-dimensional wave-guiding theory 
is still useful as a guide to the behavior of a complete theory. In the full three-dimensional case, the area term A is 
simply omitted, as A = 1. 

The resulting generalized Maxwell-Bloch equations are: 

- c^V^] A = c^V X P(x) 
[d^ + nl{k)] = g,(x)V X A + V,; (a,(x)V,;p,(x)) . (6.3) 

In order to quantize the system, we introduce the canonical momenta , n(x) and tt^{x), where: 

n(x) = M(x) = B(x) , 

7ry(x) = — p^(x) oc TOV . (6.4) 

£o5f(x) 

The quantization rules are the usual ones obtained on replacing Poisson brackets with operator commutators. 
Scaling by A is introduced so that the field units are independent of waveguide dimension, and the delta-functions are 
n-dimensional. All position-type operators of the same type of variable must commute amongst themselves, as do all 
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momentum-type operators. The commutators involving cross-terms between the position and momentum operators, 
and between the different types of variable, are: 



(I) 


■a(x),B,(x')" 


= ihVx X ^,:^(x 


(II) [Pw(x) , TTjv (x')] = ih6ijS„^'S{x - 


(III) 


A(x), Pjv(x') 


= 0, 


(IV) 


Bi{x), 7rjv(x') 


= 0, 


(V) 


A(x) , Tfj„(x') 


= , 


(VI) 


-Bi(x) , ftv(x') 


= . 



(6.5) 



Here we have introduced the usual notation of a transverse delta function defined so that: 



1 



(27r)^ 



d"k(5T^(k)e' 



ik-x 



(6.6) 



where i5,^(k) = (Sij — fcjfcj/|k|^) is the transverse projector in momentum space. In three dimensions, the first 
commutator can also be written in the more familiar form of: 



(I) 



Ei{x) , A(x') = X S^j{x - x')/£o 



(6.7) 



The final quantum Hamiltonian density is (using normal ordering): 



H/A = ^D2(x) + J-B2(x) - 1d(x) . P(x) 



+ 



1 



N 



2eo5'i^(x) 



[eg5^(x)7f2(x) + 172p2(x) + a.(x)[V.p.(x)]2] 



(6.8) 



i/=i 



These equations hold for an arbitrary spatial distribution p(x) of the continuum of polarizable atoms. 



VII. THREE-DIMENSIONAL QUANTUM MODE OPERATORS 

In this section, we specialize to the case of a continuum model with uniform couplings and velocities, as in the 
one-dimensional case, in order to find the dispersion relations for plane-waves. As in the one- dimensional case as well, 
we introduce a mode structure by finding the eigenmodes to Maxwell's equations in the form: 



A(t,x) \ / A(t^,k) 
Pj.(i,x) ) \ p^(w,k) 



Defining k = |k|, it follows that these satisfy the Maxwell-Bloch equations in the form: 

(-w^ -I- c^fc^) A = ic^k X ^ , 

{-Lo'^ + r!2(fc))p, = ig.li X A . (7.2) 

Our model now includes simple phonon/exciton dispersion. These effects cause the Fourier-domain equations for 
the polarization to be modified, so that fi^ is now replaced by a momentum dependent function fi^(A;) =Vt1 + k^a„, 
where we define k = |k| in this section. The phonon/exciton dispersion can be generalized to more complicated cases 
if desired, with higher order k-dependences - although in many cases only the relatively long wavelength (small k) 
phonons are relevant to optical properties. In terms of the coupled equations given above, only the transverse part of 
the polarization couples to the photons, hence we can define: 

pi{k) = ~6tj{k)pj{k) . (7.3) 
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Solving for P-"- by summing over the transverse polarizations of all the oscillators, we find: 



(7.4) 



The eigenvalues of the longitudinal part of the polarization field are not changed by the coupling, while eliminating 
the transverse polarization field P""" leads to the eigenvalue equation: 



n'^(uj) 



9u 



(7.5) 



For any wave-number k = |k|, we find a band-structure in which there are typically (2(iV + 1)) positive roots 
uiocr{k), ..UNaik) to the transverse eigenvalue equations for each k value, with a = 1,2, and a;^i(fc) = w^2(fc) due 
to the isotropy of our model. Each mode therefore has two orthogonal polarization unit vectors Uo^(k), such that 
k • u<T(k) = 0. Similarly, there are (N) longitudinal roots, which are labeled a = 0, and are unchanged by the 
long wavelength electromagnetic couplings. In the case where phonon/exciton dispersion is omitted, the transverse 
dispersion relation is precisely the same as in the one-dimensional model. 

As before, the solution is unique for any given modal frequency, but has forbidden regions which indicate a resonance, 
or absorption band. Typical dispersion relations for this model also demonstrate the existence of transmission and 
absorption bands - but with additional structure in each branch, whose origin is in the phonon (or exciton) dispersion. 
Prom now on, we use the notation w^cr(^) with a = 1,2, to indicate a solution to the full transverse equations. The 
notation uJ^o{k) ~ f^p(fc) indicates the longitudinal eigenvalue, which of course is undefined for n = 0, in the absence 
of longitudinal photons. 



A. Three-dimensional commutators 



Having derived the modal solutions, we now wish to expand the fields of the theory in terms of annihilation and 
creation operators in the three-dimensional model. These have the function of diagonalizing the Hamiltonian, which 
we anticipate will have the final structure of: 



N 2 



H 



U = (T = 0' 



(7.6) 



Here the lower limit notation of ct = 0' is defined to exclude the unphysical combination of /U = 0, cr = 0. As before, 
the modal frequency ujf^a{k) is the inverse of fc(il), for the /i-th transverse branch. This expansion requires that we 
define mode operators a^a in the fi-th branch of the dispersion relation so that: 



N 2 



A(i,x) = E E / d"k [A^<,(k)a^<,(k)e'''— ^--('^^^ + h.c. 

^=0 cr=l 



(7.7) 



The summation here is just over the transverse branches in the dispersion relation. These combined transverse field- 
polarization excitations can be termed polaritons, and we notice here that there can be longitudinal excitations that 
propagate, as is usual in solid-state physics. The commutation properties of the annihilation and creation operators 
(k) are chosen so that: 



a^,(k),aj^,,,(k')] = d^^,6^^,6{k - k') . 
Similarly, the momentum field can be expanded as: 

N 2 

n(t,x) = E E / d"k [n^,(k)a^,(k)e'''-^-''^''<'('=)* + h.c. 

^=0 cr=l 

Thus, at equal times: 

N 2 

A,(X) , 7f,(x')] = »M,:^(X - X')M = E E / '^"'^ 



fj,=Q a=l ' 



(7.E 



(7.9) 



(7.10) 
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This implies that, in order to preserve commutation relations, we have similar results to the one-dimensional case - 



N 2 

^^A.,.(k)n*,,(k) = ^j^s^A^) ■ (7.11) 

p=0 cr=l ^ ' 

For free fields, it is clear that 11 — fiA. Hence, we can also write n^o-Ck) in the form of: 

n^,(k) = -iLo^,{k)^IA^,{■k) . (7.12) 

The equation for the expansion coefficients A^cr(k), is therefore: 

N 2 ^ 

^^a,,.(fc)A,,.(k)A,,.(k) = / ,^(k) . (7.13) 

fj,=Ocr=l ' 

Next, we choose A^^(\i) to be real, and as before, we can always choose an (unknown) function Vf^a{k) so that: 



A^^(k) = u^(k) 



2Ak{2n) 



1/2 



(7.14) 



We will show later that v^„{k) must be interpreted as the electromagnetic component of group velocity, with non- 
electromagnetic phonon/exciton dispersion explicitly excluded. This is not completely identical to either the earlier 
narrow-band analysis of this problem, or the simple one-dimensional theory. The difference can be attributed to 
the difference in the Hamiltonian energy that is introduced when the polarization fields are dispersive. 

In order to demonstrate this, we first recall the standard identity that: 

2 

5^«,,(k)«,,(k) = ^^(k) . (7.15) 

(T = l 

Next, substituting the above expression for Apcr(k) into the equation for the field and mode commutators leads to: 



N 



(I) E^ = i- (7-16) 



As explained above, we have to determine a function v^a{k) which satisfies this condition, and we intend to 
demonstrate that the choice of v^a(k) as the (isotropic) electromagnetic component of group- velocity of the relevant 
polariton branch is sufficient - no other correction factors are needed. At this point, we notice an important fact; apart 
from the change in the resonance frequencies Oy(fc) due to their k-dependence, the above summation over ujfj^^{k) is 
identical in analytic form to our one-dimensional expression. However, the k-dependence plays no role in obtaining 
the Cauchy theorem summation results, provided we define v^„{k) to be the same function of ^^{k) and Lij^(j{k) as 
before. In other words, the k-dependence of the phonon-exciton dispersion relation simply renormalizes the effective 
resonance frequency at each k-value, in the above summation over the roots of the dispersion relation. Since this 
renormalization is the same for each eigenvalue, the summation can be carried out using identical techniques to those 
used previously. Thus, we define, for <j — 1^2: 



It should be noted that this v^^[k) is not the same as Vf^crik) in the one-dimensional case, although it has the same 
functional form. The difference is that it is a now a function of f^^(fc) and u!i^„{k), which include phonon/exciton 
dispersion effects. However, if we simply differentiate the slope of the function uj^„{k), we do not obtain as defined 
here - there are additional terms involving dfl^{k)/dk. For this reason, we refer to Wpo-(fc) as the electromagnetic 
component of group velocity, which excludes additional transport terms. 

It is clear that the mode function expansion of p,/, tTi, is also needed, for a complete demonstration of consistency. 
Using Maxwell's equations, if we expand Pi, as: 



18 



N 2 
lj,=0 cr=0 



p;:,(k)a^<,(k)e*^-^--('=)* + /i.c. 



(7.18) 



then it follows that the expansion coefficient for the transverse component of p^^ in the fi-th frequency band must 
be: 



PL(k) = 



ig^k X Apcr(k) 



Similarly, if the canonical momentum for the atomic polarization field is: 

N 2 

' (i"k 7r;;^(k)a^^(k)e''"''-''^''<'('=)* + /i.c. 



^=0 a=0 



then the corresponding transverse expansion coefficient is: 



<.(k) 



(7.19) 



(7.20) 



(7.21) 



For these operators to have the correct equal-time commutators, the diff'erent oscillator position operators Pi, must 
commute amongst themselves at equal times, as must the different momentum operators tTj,. This is trivial from the 
form of the mode operator expansion. However, the commutation relations (II) between the position and momentum 
operators are non-trivial. The relevant commutation conditions are therefore: 



N 



\pi4x),nj^'{x')]=inS^^'SijS{x-x')/A = / '^"'^ 

fl=0 cr=0' 

This in turn implies that: 



PL.(k)7r*;::(k)e^'^(---')-/i.c. 



(7.22) 



EE^'r..(k)7r;;;(k) 

IJ,=0 cr=0' 



(7.23) 



Defining 51 -{k) = kikj/k' — ^„ 



2A(27r)» 

^ — Sii — (^^^(k), and expanding the coefficients, gives two new equations. The transverse 



case is: 



(11) 



N 

E 



c^k'^v^„{k)gi, 



.{k){ujUk)-ni{k)){^l,{k)-nl,) 



(7.24) 



while the longitudinal equation is easily solved on defining Uo(k) = k/fc: 

1/2 



Pjlo(k) = uo(k)(5^^ 



2^(27r)"f2^(/c) 



(7.25) 



Finally, to ensure that there are correct field-atom commutators, we must satisfy the commutators III- VI. For these 

cross-terms between the oscillators and field variables, we find that conditions (III) and (IV), involving commutators 
between the field and the particle position (or the field momentum and particle momentum) are automatically satisfied. 
This occurs for the same reason that commutators like [A(x), A(x')] or [7r(x), tt (x')] must equal zero in our mode 

expansion. In all these cases involving pairs of canonical position-type operators or pairs of momentum- type operators, 
the commutator reduces to an odd integral over fc, which vanishes when integrated over all positive and negative k- 
values. 

This leaves the requirements (V) and (VI), which are that A and tt^ must commute at equal times , as well as tt 
and Yv These two requirements both imply the same restriction on the expansion coefficients, and hence on Vfj,a{k), 
which is that for all k and u we must have the condition: 



kv^a{k) 



N 



. 



(7.26) 
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Despite the complex nature of each of these conditions - which involve sums over all the roots of the dispersion 
equation, and must be satisfied for all the resonant frequencies w^o-, as well as all momenta k - we will show that 
each of these sums can be analytically evaluated without requiring an algebraic solution for the roots, just as before. 
For all of the commutation relation identities it is preferable to use techniques from complex function theory, which 
transform the sums over roots of the dispersion relation to complex contour integrals of related mcromorphic functions. 
However, the transverse dispersion relations considered here have an identical analytic structure for a fixed k-value, 
with those in in the one-dimensional case, so the previous analytic results follow without any further calculation. The 
main point here is that it is necessary for the 'group- velocity' coefficient to have the algebraic form given in Eq.(7.17) 
- which implies that it includes only part of the slope of the dispersion relation. 



VIII. HAMILTONIAN 



We now wish to show that when the Hamiltonian is expressed in terms of the operators apcr(k) and aj^^Ck), 
fi = 0, . . . , N and a = 0, 1, 2, where the combination /i = cr = is omitted, it is of diagonal form. In order to prove 
this we shall study further the classical modes, in particular their orthogonality and normalizat ion properties. Once 
this has been done, we shall be able to show that the Hamiltonian takes the form given in Eq. (7.6). 

Let us begin by restating the classical mode equations, Eqs. (10), in matrix form. Define the 3N + 3 component 
vector A by 



A 



(8.1) 



or Ao = A and As = pj.^ for s > 1, and the 3(iV + 1) x 3(7V + 1) matrix M(k) by 



Af(k) 



/ k^c^h -ic^K -ic^K . . . \ 
-ig^.K nl(k)h 



-ig,,K KAk)h 



V 



(8.2) 



Here is the 3x3 identity matrix and K{\<i) is the anti-hcrmitian matrix given by 

-fca fca 
K=\ k:i -ki 
-k2 ki 

which has the action on an arbitrary vector A 

KA = 'kxA. 

The equations for the modes can now be expressed as (for each value of k) 

MA = LU^X. 



(8.3) 



(8.4) 



(8.5) 



The matrix M is not hermitian, but if it is multiplied by the positive, diagonal, 3(A'^ + 1) x 3(iV -|- 1) matrix G, 

...\ 



G 



then the combination GM is hermitian. Therefore, if 




























(8.6) 



MA<^' = c.?A<^' 



u.rx''\ 



(8.7) 
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then 



(X"\GMX'") =ojl{t''\GY') (8.8) 

= (GAfA^''|A*'^) =cc;2(A^'Va^'V (8-9) 

This imphes that if uif ^ then 

(A^'Va^'^ = 0, (8.10) 

and we have part of the desired orthogonahty relation. 

In order to learn more we must examine the 3A^ + 3 eigenvectors in more detail. Define the projection operator, 
which projects each component of A onto its longitudinal component, 



P(k) 



/|k)(k| .. 

|k)(k| ... I , (8.11) 



where |k)(k| is the projection onto the vector k. A short calculation shows that [P(k),Af] = which implies that 
the eigenvectors A can be taken to lie in either the subspace projected out by P(k) (longitudinal modes), or in the 
orthogonal subspace (transverse modes). The longitudinal modes can be found by taking the inner product of Eqs. 



(7.2) with k giving 

{k^c^ -uj^)ii- A = (nl{k) - uj^)k ■ ^ 0. (8.12) 

There are N physical solutions to these equations given by A = 0, p^^ oc k, and Pi, = ior v Vf^^ with eigenvalue 
Vil (k) for — 1, . . . N. There is also one unphysical solution (it violates the gauge condition) given by A oc k and 
all of the Pi, being equal to zero. 

We are now left with 2iV+2 transverse solutions. We first note that each transverse eigenvalue is two-fold degenerate. 
This follows from the fact that if MX — uj'^X, then M{K3nX) = uj'^{K3nX), where 




i^sAT - 1^ " ^ ■ • • j , (8.13) 

which can be verified by noting that if A and Pi, satisfy Eqs. j^T^), so do k x A and k x Pi,. We choose the two 
eigenvectors A'^ and A^ , which correspond to the eigenvalue w^, to be orthogonal in the sense that 

{T^\gT^)=0. (8.14) 



This, along with Eq. ( S.IO ) implies that 



(A^^^Va*"''''^)=5.m'^-', (8.15) 

which is our final orthonormality relation. Here, due to isotropy, the mode frequency tj^cr does not depend on the 

polarization index, ct, for the transverse modes. 

We now express the fields in terms of the eigenvectors, substitute them into the Hamiltonian density in Eq. ( |6.^ ), 
and integrate over the relevant n-dimensional volume. In particular, we have that 

A^,''"^ = A^, A(^-) - p;:, (8.16) 
n^, = -^^iw^A,,, = ^^Ai^"). (8.17) 

Adding the requirement that A^^^V)* = A^^^^k) (it can be shown that M(-k)A*^''' (k)* = ujl{k)>!''"'\k)* which 



implies that A^'^'^^(k)* is in the two-dimensional subspace spanned by a'"^'^''(— k) for cr = 1, 2) and utilizing Eq. ( ^.15 ) 
we find that the terms of the form a^o-(k)a^'o-' (— k) vanish giving for the transverse modes 
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N 2 

n 1 'J 



p=0 cr=l 



In order to show that the Hamiltonian assumes the form given in Eq. (7.6) we need to prove that 

2(2^)'V^(A^''"\k)|GA^''"\k)}w^(A:) = Ti. 
Making use of Eqs. ( 7.14 ) and ( [7.19 ) this condition becomes 



1+E 



1. 



(8.18) 



(8.19) 



(8.20) 



This agrees precisely with Eq. (7.17) and is true even including phonon dispersion (a^ =/= 0). However, as pointed 
out earlier, when there is phonon dispersion we cannot interpret u^o- as the total group velocity - it only includes an 
electromagnetic contribution, i.e., it is no longer equal to dto/dk. 

In order to complete the diagonalization of the Hamiltonian we must consider the longitudinal modes. The fields 
A and 11 have no longitudinal components so that this part of the diagonalization procedure involves only the fields 
Pj, and TTy. We find that 



OA ^ r (k) 
Hiong = {2n)'— E / d'k^^p^^,{kr ■ p;io(k)4o«.o(k) 



.21) 



where we have made use of the fact that pj;^g(k) = (5^j,p^Q(k), and we note that ujf^o{k) = ri^(fc). The Hamiltonian 
assumes the expected form. 



JV 

^^long = X! / d^khu!i_,o{k)al^Qaf^o{'k), 



.22) 



when explicit expressions for the vectors, p^Q(k) from Eq. ( 7.25| ) are used. 

The final Hamiltonian, which is the sum of -fftrans and -ffiong, has (3iV + 2) mode operators for each value of k, and 
can be written in the form 



N 2 

^ = E E / ^cu^.{k)aU^)a^A^)rk 

U=0 (7=0' 



(8.23) 



Here the lower limit cr = 0' excludes the combination of ^ — and cr = 0, which would imply a longitudinal polariton. 
Also, there is a requirement of having (iV + 1) distinct roots for this form to be valid. The corresponding field operators 
(in the full three-dimensional case) have the expansions: 



D 



B 



^ E 

M,(T = 1,2 

^ E 

M,(T=1,2 

- E 

fj.,a=l,2 



d^k 



nkv^'"'{k)eiu^{k)) 



-.1/2 



hkvf;''{k) 



ecr(k)a^^(k)e 



ik-x 



he 



1/2 



d^'k 



47re(u;^(fc)) 
h^ikv;^"{ky^^^ 



ea{k)ai,^{k)e"^'' + hc 



47r 



Uaik)af,aik)e'^''' + he 



.24) 



Here we have introduced the electric field mode e(k) = k x u(k)/|k| to simplify the expansion. We note, as in the 
one-dimensional case, that the transverse field expansion for the electric field is simply derived from the displacement 
field by using the frequency-dependent permittivity. The main feature introduced by the dispersion is the replacement 
of a frequency term u, that would normally appear in the expansion coefficients, by a new term with the same units, 
but equal to kv^'"{k) instead. We do not give the expansion for the longitudinal part of the electric field here explicitly, 
except to point out that it is equal to — P" /eq. 
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IX. SUMMARY 



A simple theory of a one-dimensional, dispersive waveguide was introduced, including a polarizable model of the 

medium with N discrete localized resonances. This can be thought of as a limiting case of an ideal insulator, in which 
the polarization field is due to localized electrons at each atomic location. The theory is exactly equivalent to the 
usual classical theory of a dispersive dielectric medium, in the sense that it results in the Sellmeir equations for the 
refractive index. These are well-known to lead to an excellent fit to the classical dispersion properties of transparent 
media, and have the usual causality requirements automatically satisfied. The theory was quantized and a set of 
+ 1 mode operators introduced, for the polaritons in each branch of the dispersion relation, provided there were 
N + 1 distinct, positive roots. In this case, the mode-expansion has a universal and simple form, only depending on 
the group velocity. 

This model is necessarily causal, and implements the causality requirements through band-gaps, rather than isolated 
poles. It does omit many important correction factors that occur in practice. In particular, our mode expansion 
neglects scattering off inhomogeneities. For this reason, transmission inside the transmission band is essentially 
lossless. It also omits nonlinearities due to phonon-phonon, photon-photon, and photon-phonon interactions, which 
are responsible for additional non-electromagnetic damping of the polaritons. However, these effects can certainly be 
added to the Hamiltonian once a mode expansion is established. 

Next, a quantum theory of an isotropic n-dimensional dispersive waveguide was introduced, with n = 1,2 or 3. 
Without any additional phonon/cxciton dispersion, the theory is exactly equivalent to the classical Drude-Lorentz 
theory of a dispersive dielectric medium. The complete n-dimensional theory was quantized, and a set of (3A^ -|- 2) 
mode operators introduced, for each branch of the dispersion relation; again, with the restriction of distinct, positive 
roots. As in the one-dimensional case, the mode expansion depends on the permittivity and the electromagnetic group 
velocity, in the case of transverse polaritons. However, the group velocity factor in this case is modified to include 
only the electromagnetic component of the group velocity. 
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